(navigation image)
Home American Libraries | Canadian Libraries | Universal Library | Community Texts | Project Gutenberg | Biodiversity Heritage Library | Children's Library | Additional Collections
Search: Advanced Search
Anonymous User (login or join us)
Upload
See other formats

Full text of "COURS_CMN_3EME_ANNEE"

Chapitre 3 
Systemes d'equations algebriques 



Mathematiques-lnformatique, Sciences de I'Homme et de la 

Societe 

Universite Moulay Ismail - Meknes 

1 1 decembre 201 1 



Q Introduction 

Q Elimination de Gauss 

Q Decomposition LU 

Q Effets de I'arithmetique flottante 

Q Conditionnementd'une matrice 

Q Methodes iteratives pour la solution d'un systeme lineaire 



Introduction Introduction 



On a resolu une equation de la forme 

f(x) = ou x e R et f(x) e R (i.e. f : R -k R) 



On complique encore les choses : deux grandes families de problemes 
possibles : 

9 les systemes d'equations lineaires de n equations et n inconnues : 

Ax = b 

ou A est une matrice carree inversible et b un vecteur 

(F : R n ->: R n avec F(x) = Ax - b et on cherche F(x) = 0) 

9 les systemes non lineaires : F(x) = ou F : R n ->■: M n est une 
application a n variables : a? equations non-lineaires dependant de 
n variables. 



Introduction Introduction 



Beaucoup d'applications dans le domaine de I'ingenierie exigent de 
resoudre des systemes d'equations qui comportent souvent un grand 
nombre d'inconnues. On rencontre ces applications, par exemple, 
dans les domaines suivants : 

» la mecanique des fluides : ecoulement de I'air autour d'un avion, 
ecoulement dans une turbine hydraulique, diffusion dans les 
arteres d'un medicament enrobant des stents, etc. 

9 la mecanique des solides : deformation d'un materiau soumis a 
des contraintes exterieures, analyse des structures complexes, 
etc. 

9 Optimisation de la distribution des frequences sur des antennes, 
filtrage de signaux, etc. 



Introduction Introduction 



Comportement hygro-mecanique de lames de parquet en services 
(periode de 42 jours, humidite relative augmentee de 20%). 




Ce calcul exige de resoudre un premier systeme lineaire de 35890 equations 
pour le mouvement de I'humidite puis un systeme de 107670 equations pour 
la deformation. Ces deux resolutions sont faites au total 129 fois pour decrire 
le comportement a divers moment dans le temps. 



Un systeme d'equations lineaires est un systeme de n equations 
lineaires a n inconnues 

anxi + a 12 x 2 + a 13 x 3 + ... + a 1n x, 7 = *>1 
a 2 i x, + a 22 x 2 + a 23*3 + • • • + a 2n Xn = b 2 



a n1 Xt + a n2 x 2 + a n3 x 3 + ... + a nn x n = b n 
ou a/y et by sont connus pour /,/ = 1 , ..., n. 



Exemple : Trouver x^ , x 2 , x 3 solution de 



2x-i + 3x 2 - x 3 = 5 
4xi + 4x 2 - 3x 3 = 3 
-2xt + 3x 2 - x 3 = 1 



Ecriture matricielle : 



Ax = b 



ou A denote la matrice carree d'ordre n des coefficients du systeme 
lineaire et b le vecteur du second membre. 



/ an 


• • a Vl 


• • a 1n \ 


a/i 


.. a,,- 


3 in 


\ a n1 


a n j 


3-nn ) 



( b, \ 



bi 



\bn J 



( *\ \ 



\x n J 



Exemple : Trouver xg! 3 solution de 



2 3 
4 4 
-2 3 



-1 
-3 
-1 



*1 

^3 



5 

3 
1 



Rappel d'algebre 



Trois enonces equivalents : 

9 Le systeme Ax = b admet une solution unique si et seulement si 
la matrice A est de rang maximal (autant d'inconnues que 
d'equations et pas d'equation qui se repete). 

9 Le systeme Ax = b admet une solution unique si et seulement si 
det/\ / 

» Le systeme Ax = b admet une solution unique si et seulement si 
A est inversible. 

Dans ce cas la solution est fournie par 

A~ 1 b 

detA 



9 X 

9 La formule de Cramer 



det4 



ou Aj est la matrice A dont la colonne i est remplacee par b. 



Remarque 



a La plupart du temps on traite des cas ou la matrice est inversible 
(matrice non-singuliere). Ce qui assure I'existence d'une solution. 

9 Sauf dans de rare cas la formule de Cramer peut etre considerer 
comme purementtheorique (sans applications). 

9 Sauf dans des cas particuliers il est plus rapide et simple de 
resoudre le systeme lineaire que d'inverser la matrice (meme "a la 
main"). 



Comment resoudre efficacement des systemes lineaires? 



Methode des substitutions successives 



La methode la plus simpliste (et classique) consiste a eliminer des 
inconnues dans les equations par substitution successive : 



2 3 

3 4 



De la premiere equation on a 



*1 



x 2 



8-3x 2 



8 
11 



Et en substituant dans la deuxieme (eliminant x^ona 

1 1 => x 2 = 2 



3(^-^)+4x 2 



Elimination de Gauss Matrices diagonales 



Pour une matrice quelconque cette approche est simple pour un 
humain mais difficile a traduire en algorithme. On doit abandonner 
cette approche du moins sous cette forme. 

Cette approche est plus simple a decrire pour certaines structures de 
matrices. 



Matrice diagonale 



Facile lorsque le systeme lineaire a une matrice diagonale : 

xi \ / £>i 

bn 



Ax 



I an 





322 


. 
. 


. 


U 




. 


. 



^>aiiXi = bi i=-\,..., 



b n 

alors X; = —. Solution unique ssi a„ / V/ (6e\A = f\ a„ / 0) 

a n /=1 



Definition 3.1 



Une matrice est dite triangulaire inferieure (ou superieure) si toutes 
ces entrees a,y sont nulles pour / < j (i > j resp.). 

Triangulaire inferieure = "nulle au dessus de la diagonale". 
Triangulaire superieure = "nulle au dessous de la diagonale". 



/ an 

321 

331 





a 2 2 

a32 






a33 



\ a n i a„2 a n3 



\ / an ai2 3i3 
a22 a23 
a 33 



V 



am \ 

a2n 
a3n 



Matrice triangulaire 



Les matrices triangulaires : resoudre sucessivement en choisissant 
I'equation avec le moins d'inconnues. 

» Matrice triangulaire inferieure : descente triangulaire (on 
commence a 1 on fini a n "en bas"). 



;-1 



x^ =— Xj 



bj - ^2 a ik x k 



/<=1 



i = 2,...,n 



a Matrice triangulaire superieure : remontee triangulaire (on 
commence a n on fini a 1 "en haut"). 



b n 

X n - — Xj 



bj- Y^ a ik x k 

k=i+^ 



i = n - 1, n - 2, ...,1 



, , , ■_,..., 



Math-lnfo-ScHS (ENSAM) 



Unicite dans le cas des matrices triangulaires 



La remontee et la descente triangulaire n'ont de sens qu'a la condition 
que les termes diagonaux soient tous non nuls. 

Que se passe-t'il si cela n'est pas le cas ? 

Puisque le determinant d'une matrice triangulaire est le produit des 
termes de la diagonale, alors 



det/4 = Y[ an ^ <£> a,/ ^ V/ 



/=i 



Done I'unicite de la solution impose que les termes diagonaux soient 
tous non nuls. Dans le cas contraire le systeme n'a pas une solution 
unique et A n'est pas inversible. 



Exemple de resolution 



Exemple de remontee triangulaire 



Soit a resoudre le systeme 



3 1 3 \ / X! 

2 2 x 2 

1 / V x 3 




30 
15 



x 3 = — = 15 x 2 

^33 



b 2 - V a 2k x k 

k=3 _ 30 - 2x 3 

322 2 



fr) -J^auXfc 



*1 



ft=2 

an 



- 1 x 2 - 3x 3 



x 3 = -15 



L'algorithme de remontee (ou de descente) est efficace. On sait que 
Ton peut transformer un systeme quelconque en systeme triangulaire. 

On voudrait 

9 une methode de passage d'un systeme de quelconque vers un 
systeme triangulaire tout en conservant la solution du systeme 
original. 

9 utiliser la remontee ou descente triangulaire. 

Pour cela on va revoir les operations elementaires sur les systemes 
lineaires. On produira ensuite deux methodes directes : la methode 
de Gauss et la methode de la decomposition LU. 



Definition 3.2 



Une methode de resolution d'un systeme lineaire est dite directe si la 
solution du systeme peut etre obtenue en un nombre fini et 
predetermine d'operations. 



Remarque sur les methodes directes 



Une methode directes est essentiellement une methode basee 
sur la remontee (descente). Pour de telle methode on peut faire 
un estime a priori du nombre d'operations en virgule flottante et 
du temps de calcul (une remontee ou une descente exige « rf 
operations elementaires). La plupart du temps la limitation dans 
I'utilisation de methodes directes vient de I'espace memoire requis 
pour la resolution. 

A I'opposee de ces methodes on trouve les methodes iteratives. 
Ces methodes ne garantissent pas la convergence et ne peuvent 
garantir un nombre fixe d'operations pour I'obtention d'une 
solution. Cependant la plupart des methodes iteratives exigent 
moins d'espace memoire. 



Transformation = produit par une matrice 



On veut garder la linearite du systeme alors on transforme avec 
les operations de bases (+, -, x, -=-) avec des scalaires sur et 
entre les lignes du systeme. 

Une transformation de ce type est representable par une matrice 
multipliant la matrice A et le vecteur b du systeme : 

WAx = Wb 



9 Quels conditions sur W? Si on a la solution (unique) de 
WAx = Wb alors on veut que cette solution soit la solution de 
Ax = b. Alors il faut que la matrice W soit inversible ! 

9 On s'assurera de la validite de notre transformation en s'assurant 
que son determinant est non nul (ou quel est inversible). 



Exemple 



Soit 



Ax 




WA = 



EH— 






Wb = 



la solution x = (3 2 1 ) r correspond a la solution du systeme original. 



W 



I 



— I 

1 


— 1 

-1 












WA 




Wb 



Le systeme n'est plus valide : pas unicite done pas correspondance 
entre les deux problemes ! 



Math-lnfo-ScHS (ENSAM) 



Operations elementaires sur les lignes 



Pour reduire un systeme lineaire a la forme triangulaire on utilise 3 
operations elementaires sur les lignes 

9 Multiplication d'une ligne /., par une scalaire A non nul : 

XL, — > L, 



9 Addition de deux lignes : 



9 permuter deux lignes : 



L, + XL] — > L, 



Li ^ Lj 



On sait qu'en multipliant (a gauche) le systeme original par une 
matrice inversible la solution sera inchangee. On va montrer que les 3 
operations proposees correspondent a des produits par une matrice 
inversible, nous assurant ainsi de conserver la solution originale. 



Operations sur les lignes : matrice f 



A I'operation multiplication de la ligne / par A non nul : A/., — y L, 
correspond le produit de Ax = b avec la matrice elementaire M : 



m- n 



1 j= 1,...,/- 1,1 + 1,...,n mn = A mj k = sinon 



On sait que det/W = A done la matrice est inversible si A / 



M 



( 1 o 

o '■■ 

... 


\o ... 





A 



o \ 




'1' / 



M- 1 = 



/ 1 o 

o '■■ 

... 


\o ... 



\ 






Operations sur les lignes : matrice P 



A I'operation intervertir deux lignes du systeme : L t < — > Lj correspond 

le produit de Ax = b avec la matrice elementaire P (matrice de 

permutation) obtenue en permutant les lignes / ety de la matrice 

identite 

/ 1 \ 

"■ 

P= j^ .. 
'-► 

\o 1 J 

Puisque PP = I \a matrice P est inversible et son inverse P~ 1 = P. On 
a aussi que detP = -1 














1 





1 








Operations sur les lignes : matrice " 



A I'operation addition de deux lignes L, + XL] — ► /_, correspond le 
produit de Ax = b avec la matrice elementaire T 



tkk = \ k = 1 , ..., n tjj = A tfd = sinon 

On sait que detT = 1 alors 7 est inversible (il suffit d'enlever ce que 
Ton a ajouter) 



/ 1 

o '■■ 

... 


\0 ... 





1 

A 



\ 






/ 1 o 

o '■■ 

o ... 
o 

\o ... 



... 

1 ... 

-A 1 



\ 






1 / 



La methode d'elimination de Gauss consiste a utiliser les operations 
elementaires afin de reduire le systeme sous la forme triangulaire 
superieure. 

On distingue deux cas : 

9 Sans pivotement : on ne permet pas la permutation de deux 
lignes (la transformation P est interdite) ; dans certain cas on ne 
pourra pas obtenir de matrice triangulaire sup. ! 

9 Avec pivotement : on permet la permutation des lignes ; dans ce 
cas on obtiendra toujours une matrice triangulaire sup. si la 
matrice est inversible. 

Dans la pratique on ne construit pas les matrices (de types M,T et P) 
necessaires a la transformation. 

Cependant la methode de Gauss est valide parce qu'elle consiste a 
multiplier le systeme de depart par une suite de matrices inversibles. 
On verra plus tard comment ces matrices nous seront utiles. 



Dans I'application de la methode de Gauss on precede par elimination 
des termes sous la diagonale en utilisant la notion de matrice 
augmentee. 



Definition 



La matrice augmentee du systeme Ax = b est la matrice de dimension 
n x n+1 obtenu en ajoutant le membre de droite b a la matrice A : 



" an . 


• a Y] 


a-\ n 


&1 " 


a/1 . 


■ an 


a i n 


bi 


. 3/71 • 


ani 


3nn 


bn _ 



Cette notation est tres utile puisque les operations elementaires 
doivent etre faites sur la matrice A et sur b simultanement. 



9 Un systeme Ax = b a une solution unique ssi 

* 6etA ^ 

• A est inversible (non singuliere) 

» # equations = # inconnues et pas dequations repetees 

9 Les systemes diagonaux et trinagulaires (inf. et sup.) peuvent etre 
resolus par "remontee" ou "descente" triangulaire. 

• Transformation d'un systeme quelconque en systeme 
triangulaire : on veut passer de Ax = b vers 

Ax = b 



avec 

a 7l triangulaire 



» la solution des deux systeme doit etre la meme 

Ax* = b ■& Ax* = b 



On a montrer que les bonnes transformations sont des produits 
par des matrices inversibles et que les 3 operations de bases : 

* multiplication d'une ligne par un scalaire mene a une matrice M 

inversible 
o addition de deux lignes mene a une matrice T inversible 
» permutation de deux lignes mene a une matrice P inversible 

Avec ces trois operations on peut transformer un systeme en 
systeme triangulaire. 

Methode de Gauss : utilisation de combinaisons de M, T et P 
pour obtenir une matrice triangulaire superieure. Deux choix : 

» Pas de permutations permises : peux mener a des matrices 
non-triangulaires 

» avec permutation : on peut toujours obtenir une matrice traingulaire. 

Dans la pratique cette methode ne tient pas compte des matrices 
de transformation. Elles ne sont la que pour demontrer que la 
methode de Gauss transforme de la bonne maniere le systeme de 
depart. 



Math-lnfo-ScHS (ENSAM) 



Exemple 3.8 



On applique la methode de Gauss sur la matrice augmentee 



"212 


10 1 


6 4 


26 


8 5 1 


35 



Avec une permutation : Exemple 3.9 



On applique la methode de Gauss sur la matrice augmentee 



1 


1 


2 


1 


2 


2 


2 


5 


3 


4 


1 


3 


3 


3 


-2 


1 


1 


4 


5 


-2 



Et si je veux changer b ? 



Trois options : 

Q Je recommence le processus avec la nouvelle valeur de b car je 
dois transformer A et b simultanement dans la methode de Gauss 

O j'attend d'avoir TOUT les b que je veux et je fais Gauss sur la 
matrice augmentee de tout les b 

Q J'essaie de garder I'information ayant servi a la transsformation de 
la matrice lors de la premiere rsolution pour resoudre avec les 
autres b 

Revenons a nos exemples pour voir ce qu'on peut faire en tenant 
compte des matrices T,M et P cette fois 



Exemple 3.8 revisite 



On applique la methode de Gauss sur la matrice augmentee 

A 



"212 


10 1 


6 4 


26 


8 5 1 


35 



Mais on construit les 3 matrices 7 et on a 

A = (T-'T-'T-')U = LU 



Avec une permutation : Exemple 3.9 revisite 



On applique la methode de Gauss sur la matrice augmentee 



1 


1 


2 


1 


2 1 


2 


2 


5 


3 


4 


1 


3 


3 


3 


-2 


1 


1 


4 


5 


-2 



Dans ce cas on a des 7 et un P et on a 

PA = (P 4 7"f 1 T,T 1 T 3 " 1 P 4 T^)U = LU^A = PLU 



La decomposition LU permet de resoudre avec differents b 



Cette decomposition donne toute I'information necessaire pour 
resoudre a "repetition". Puisque L et U sont construit on a 

Ax = b o (LU)x = b^ L(Ux) = b 

Si on pose y = Ux alors la resolution de Ax = b peut se faire en deux 
temps : 

O resoudre Ly = b par une descente triangulaire 

O resoudre Ux = y par une remontee triangulaire 

Une fois le probleme resolu pour un b, en ayant garder L\e peux 
changer le b et resoudre a nouveau sans refaire de transformation du 
sy steme. 



Remarques 



9 la methode de Gauss est une decomposition LU ou on "jete" L. 
Plus la peine d'en parler (sauf pour les exercices :-) 

9 les permutations causent des problemes (empechent d'avoir un 
LU "propre") on en reparlera. 

9 Pour le moment on va faire comme si A etait permuter avant qu'on 
arrive... 



Unicite de LU 



La decomposition LU se fait grace aux transformations de Ax = b mais 
on pourrait ajouter des transformations (superflues). Par exemple on 
pourrait finir en s'assurant que les pivots (termes diagonaux non nuls) 
soit tous egaux a 1 en appliquant des transformations M. Done la 
decomposition LU n'est pas unique. 



Pour eviter les ambiguites possible (I'operateur == est mal defini) on 
voudrait assurer I'unicite de la decomposition. Pour se faire on va 
imposer une conditions supplemental a notre decomposition. 
Deux methodes "classiques" 

9 Decomposition de type Matlab : /.,, = 1 V/' (pas de division dans 
la descente) 

9 Decomposition de Crout : L/„ = 1 V/ (pas de division dans la 
remontee) 



Decomposition de Crout 



Algorithme compose de trois parties 
O factorisation LU avec U„ = 1 V/' 
O Resolution de Ly = b par descente triangulaire 
O Resolution de Ux = y par remontee triangulaire 



Comment on fait la factorisation 



On identifie termes a termes les composantes de L et U. A noter que 
c'est un systeme lineaire de n 2 equations et inconnues qui a une 
solution unique quand A est inversibles. 

Une particularity de I'approche est que Ton resout alternativement 
pour une colonne de L puis pour une ligne de U en commengant par la 
colonne 1 done on construit : 



Lft , U-\j, ..., Lj n _-\ , U n --\j, Li, 



Cette strategie nous permettra de reduire I'espace memoire lors de la 
factorisation. 



9 Factorisation 

• L;i=A- 1 i=-\,...,n,Uu = T$ i = 2,.,.,n 
» i = 2,...,n- 1 

;-i 
a L„ = A„ - J2 L ik U ki 

k=\ 
<* /=/+1,...,n 

l-ii^Aji-'^LjKUu 
/t=i 

;'-1 
A;- £ i-ftUltf 
//.. — *= 1 

n-1 

• '-nn = A nn - Y, Lnk Ukn 

/c=1 

» descente sur Ly = b 



• yi = t y>=- 

9 remonte sur Ux = y 



u, 



/ = 2,...,n 



»x n = y n x,=y,- J] I/*** ;' = n-1,...,1 

fe=;+i 



Remarque 



L'algo ne fonctionne que si /.„ / V/. 
• Si on a pas fait de permutations alors det/\ = detLdetL/ mais 



detL/ = 1 et 



deM /0« detL / o /.,-, / V/ 



• Si on a fait des permutations det/4 = detPdetLdetL/ Puisque 
detP / alors 

deM / o detL / o L, y / V/ 

Done si /4 est inversible I'algorithme va bien. 



Remarque 



D'un point de vue pratique une fois la factorisation LU faite on a plus 
de raison de conserver la matrice A. En fait notre procedure de 
factorisation est faite pour que Ton puisse "ecraser" A au fur et a 
mesure de la factorisation. Ce qui permet de minimiser I'espace de 
stockage puisqu'on a pas simultanement A, Lei U 



Definition 3.4 



La notation compacte de LU consiste a remplacer la matrice A par la 
matrice 



Les coefficients de la diagonale de U ne sont pas garde puisqu'ils sont 
tous egaux a 1 mais on ne les neglige pas pour autant. 

Cette notation revient a stocker : 

L+{U-I) = L+U-I 
a la place de A 



Exemple 3.11 



Utilisation de la decomposition LU avec notation compacte pour 



3 
1 
2 



-1 
2 

2 



2 " 




' 12 " 


3 


b = 


11 


1 




2 



Permutations 



On gardera un vecteur de permutation (position ou numerotation) qui 
indiquera les permutations effectuees sur la matrice. A la fin de la 
factorisation on obtient L, U et le vecteur de permutation. 

Si necessaire on reconstruit la matrice P en appliquant les 
permutations sur une matrice identie. Mais c'est normalement inutile. 

Comment resoudre une fois la factorisation terminee ? 

Ax = b et PA = LU o PAx = LUx = Pbo Ly = PbeXUx = y 

On ne construit pas P mais on permutent b avec le vecteur de 
permutation puis on resoud par descente et remontee. 

Voir exemple 3.12 

ATTENTION P est en general le produit de plusieurs permutations 



P=PpPp-A..r i 



Pi o P~ 



PlPp...Pn 



Math-lnfo-ScHS (ENSAM) 



Dans le cas d'une matrice symetrique il est possible deconomiser de 
I'espace en imposant une condition supplemental a la decomposition 
LU. On impose que 

U = L T => A = LL T 

A etant symetrique on ne stocke que sa partie inferieure et sa 
diagonale. Mors on peut faire la meme chose sur la decomposition et 
ne stocke que L. 
La construction est basee sur Crout. 



Algorithme de Choleski 



• Z_n =(A u )-2 



9 /c = 2,...,A7 i = k + 1 , ..., A? 



fc-1 



• /.* 



(4-E# 



7=1 

fc-1 



-ft 



ts 



Remarques 



» Tout comme pour la decomposition LU on peut introduire une 
notation compacte et "ecraser" A lors de la creation de L 

9 La factorisation n'est pas unique : A = (-/.)(-//) = LU on rend 
la factorisation unique en imposant des pivots positifs : /.„ > 0. 



Exemple 3.13 



Utilisation de la factorisation de Choleski avec notation compacte pour 



4 
6 

2 



6 

10 
5 



2 

5 

14 



Existence de la factorisation de Choleski 



Soit 



-1 2 
2 6 



On ne peut pas faire Choleski car L u = ^/Ay\ n'est pas reels. 

En fait on aura pas existence de la factorisation si une des racines (les 
pivots) n'est pas reelles 



Definition 3.5 



Une matrice A est definie positive ssi 

Ax ■ x > Vx / 



a— 



Condition pour I'existence de la factorisation de Choleski 



On peut montrerque 

A symetrique definie positive o A a une factorisation de Choleski 



La verification de la positivite est difficile, la meilleur approche : 
appliquer la factorisation... 

Dans le cas ou on n'a pas la positivite de la matrice on applique 
simplement la factorisation LU 



9 A partir de I'elimination de Gauss on a cree la methode de 
decomposition LU 

• U est la matrice obtenue par Gauss et L est composee de 
I'information jetee par Gauss pour transformer la matrice en 
triangulaire. 

9 Le systeme se resout maintenant avec une descente (Ly = b) 
suivie d'une remontee (Ux = y). 

9 la decomposition n'est pas unique 

a factorisation Doolittle (Matlab) : /.„• = 1 

a factorisation Crout : Un = 1 $= C'est elle que I'on retient 

9 Notation compacte : on ecrase A en la remplagant par L+ U - 

9 Les permutations sont executees apres chaque construction 
d'une colonne de L et sont stockees dans un vecteur. 



PA = LU 



9 En utilisant les proprietes des matrices triangulaires et de la 
factorisation le determinant de la matrice A : 

deW = (detP r )(det/.)(detlO = (detP r )(detZ.) 

Puisque detP 7 = (-1 ) p ou p est le nombre de permutations 
effectuees on a 

n 

/=i 

9 Lors d'une factorisation le nombre d'operations x/-=- est de I'ordre 
de n 3 /3 

9 Lors d'une remontee ou descente le nombre d'operations x/-=- est 
de I'ordre de n 2 

9 On cherche a optimiser encore I'algorithme : pour certains cas ou 
la matrice est symetrique (A = A T ) on peut modifier la 
decomposition en imposant que U = L T . C'est la factorisation de 
Choleski qui entraine une reduction en espace memoire et un 
algorithme simplifie. 



La factorisation de Choleski est unique a un signe pres. Pour la 
rendre unique on impose /.„ > 

La factorisation de Choleski existe a une condition (3 versions 
equivalentes) : 

a A est symetrique et definie positive : A — A T et Ax ■ x > Vx / 0. 
• Les determinant de A et des sous-matrices principales de A sont 

tous strictement positifs. 
o Les valeurs propres de A sont toutes reelles et strictement 
positives. 

Dans le monde reel (n > 10000) la meilleure fagon de verifier 
qu'une factorisation de Choleski existe est de I'essayer. Les 3 
conditions presentees ici sont plus couteuse (numeriquement) 
que d'essayer Choleski et de ne pas reussir 

Si la factorisation de Choleski n'existe pas on fait un LU ordinaire 
ce qui implique que Ton ne profite pas de la structure particuliere 
deA. 



Exemple d'existence de la factorisation de Choleski 



Soit 

4 6 2 

6 10 5 

2 5 14 

On calcule les determinants des sous-matrices principales. 



Exemple d'existence de la factorisation de Choleski 



Soit 



-1 2 
2 6 



AlorsdeMn 



-1 < et la factorisation de Choleski n'existe pas. 



II y a une autre structure de matrice souvent rencontree et pouvant 
mener a des simplification algorithmique et d'espace memoire. II s'agit 
de la matrice tridiagonales. 



Definition 3.6 



Une matrice est tridiagonale si ses seuls termes non nuls sont situes 
sur la diagonales principale et sur les deux diagonales adjacentes (au 
dessus et au dessous de la diago. principale) 



a-11 a-12 





^21 ^22 


323 


a 32 


333 



a nn _! 







3n-ln 
&nn 



D^ 


S^ 










h 


D 2 


s 2 










k 


D 3 






S/7-1 











/ n -i 


D n 



Concernant les tridiagonales 



On stocke en memoire les vecteurs D, I, Se\ b prenant au total 
seulement 4a? - 2 reels. On fait encore plus d'economie si la 
matrice est symetrique (on ignore S dans ce cas). 

L'algorithme de decomposition est beaucoup plus simple. On 
I'obtient en appliquant la decomposition de Crout pour ce cas 
particulier. 

On retrouve cette structure de matrice dans plusieurs applications 
des methodes numeriques : 

a differences finies (approximation de la solution d'equation 

differentielle ordinaire) 
• splines cubiques (approximation par une cubique d'une fonction a 
partir d'un groupe de points) connus 



Retour sur f inverse de A 



Comment inverser A ? 



La methode usuelle (i.e. a la main) consiste a augmenter la matrice 
avec la matrice identite et a faire une elimination jusqu'a ce qu'on ait 
une "identite a gauche" 







•v~> 






A 


1 


1 


A 1 











En utilisant la decomposition LU cela revient a resoudre n fois le 

systeme Ax = b avec b = e, (les vecteurs de la base euclidienne) pour 

i=1,...,n. 

Clairement il ne taut pas utiliser cette approche pour resoudre un 

probleme : on ferait n resolution et un produit matrice vecteur. 



Combien ga coute? 



Nombre d'operations x, 



Si on veut absolumment calculer I'inverse. II faudra faire une 
factorisation LU suivi de n descente et de n remontee. Utilisant les 
ordres de grandeurs sur le nombre d'operations x/^ona 



n 


Ax = b 


A~^b 


10 


333 


1 333 


100 


333 333 


1 333 333 


1000 


333 333 333 


1 333 333 333 



Pour un systeme 1000 x 1000 on fait 10 9 operations x/-=- de plus pour 
resoudre le systeme. 



1000 est un petit chiffre dans le monde reel... 



Exemple 3.15 



On se donne une matrice A obtenu par factorisation LU avec 
permutation : 



2 1" 




1 


L = 


3 1 





3 

2 

1 







-1/3 



Avec le vecteur de permutations 



O 



3 
1 
2 



U 



10 1/3 
11/2 
1 



On va calculer I'inverse. II s'agit SURTOUT d'un exemple de resolution 
avec un vecteur de permutation 



Effets de I'arithmetiqu 



Propagation d'errt 



Exemple 3.16 1 


Soit 

A = 


1.012 -2.132 3.104" 
-2.132 4.096 -7.013 
3.104 -7.013 0.014 




On veut appliquer la factorisation LU mais avec une mantisse fixe de 4 
ch iff res. 


Cette restriction engendre des "perturbations" par rapport a la 
factorisation "humaine" (mantisse infinie). Quel est I'impact de cette 
perturbation lors de la resolution ? 



Effets de I'arithmetiqu 



Propagation d'errt 



Exemple 3.16 suite 


Avec virgule flottante a 4 chiffres 




L+U-l = 


1.012 -2.107 

-2.132 -0.3960 

3.104 -0.4730 


3.067 " 
1.197 
-8.940 




Avec virgule flottante humaine : 




L+U-l = 


1.012 -2.106719... 

-2.132 -0.395525... 

3.104 -0.473743... 


3.067193... " 

1.197757... 

-8.939140 




Quel sont les effets sur la solution de Ax = b pour b donne ? 



Effets de I'arithmetiqu 



Propagation d'errt 



Effet des operations arithmetiques en virgules flottantes 



On retourne aux considerations du chapitre 1 . Pour le resoudre un 
systeme Ax = b\\ est a prevoir que des erreurs de troncatures vont 
apparaitre. 

On verra que pour certaines matrices les effets sur la solution seront 
negligeables. Pour certaines autres matrices les variations dans les 
entrees de la matrice induiront d'importantes erreurs dans la solution 
obtenues. De telles matrices sont dites mal conditionnees 



Effets de farithmetiqu 



Propagation d'errt 



Exemple 3.17 



Soit une matrice A et une perturbation A de A 



1. 2 ' 
1.1 2 


b = 


10 ' 
10.4 


A = 


1. 2 " 
1.05 2 



La solution exacte avec A est (4, 3) et pour A la solution est (8, 1 ). 
Une "petite" variation de A produit une "grande" variation de la solution. 

Puisque l'arithmetique flottante (mantisse finie) menent a de petites 
perturbations dans la decomposition LU. II est possible que 
d'importants effets apparaissent sur la solution du systeme. 

On peut produire d'importantes erreurs sur la solution numeriques via 
la decomposition LU et le precede de descente et remontee. 



Effets de I'arithmetiqu 



Propagation d'errt 



Recettes simples pour eviter les degats 



Dependant de la matrice A on pourra avoir ou non de large erreurs. On 
verra par la suite comment detecter de telle matrice. Pour le moment 
voici deux corrections simples a appliquer lors de la resolution pour 
eviter deux des principales causes d'erreurs : 

9 division par des nombres trop petits (pivots trop petits) 

• la soustraction/addition avec des nombres ayant des ordres de 
grandeurs trop eloignees 



Effets de I'arithmetiqu 



Propagation d'errt 



Exemple 3.18 



Soit 



0.0003 3.0000 
1.0000 1.0000 



2.0001 
1.0000 



La solution est x = [1 /3, 2/3] T , en faisant la decomposition LU avec 4 
chiffres on a 



L+U-l 



0.3000 x 10" 3 
0.1000 x 10 1 



0.1000 x 10 5 
-0.9999 x 10 4 



0.6667 x 10 4 
0.6667 x 10° 



La solution par LU est "loin" de la solution reelle. Comment eviter cette 
erreur? 

On fait une division "presque par zero" a cause de An c'est lui qui 
modifie L/ 12 et L 2 2- Mors on va eviter cette division en faisant une 
permutation. 



Effets de I'arithmetique flottante Propagation d'erreurs 



Encore des permutations. 



Soit 



1.0000 1.0000 
0.0003 3.0000 



1.0000 
2.0001 



0.1000 x 10 1 0.1000 x 10 1 
0.3000 x10" 3 0.3000 x10 1 



0.3333 x 10° 
0.6667 x 10° 



en faisant la decomposition LU avec 4 chiffres on a 

L+U-l 

La solution par LU est "proche" de la solution reelle. 



Recette 1 



On peut eviter une partie des erreurs en faisant une permutation 
systematique des lignes lors de la decomposition : on s'assure d'avoir 
TOUJOURS le pivot de plus grand en valeur absolue sur la diagonale. 



Effets de farithmetiqu 



Propagation d'errt 



Exemple 3.19 



1 6 9 

2 1 2 

3 6 9 



Effets de farithmetiqu 



Propagation d'errt 



Exemple 3.20 



2.0000 100000 
1.0000 1.0000 



100000 
2.0000 



La solution est x = [1 .00002, 0.99998] T , en faisant la decomposition 
LU avec 4 chiffres on a 



L+U-l 



0.2000 x 10 1 
0.1000 x 10 1 



0.5000 x 10 5 
-0.5000 x 10 5 



0.0000 x 10° 
0.1000 x 10 1 



La solution par LU est "loin" de la solution reelle. Comment eviter cette 
erreur? 

En calculant L 2 z on fait une soustraction de chiffres d'ordre de 
grandeur tres differents (operations a eviter Chap 1) a cause de /4 12 . 
Mors on va eviter cette operations en multipliant la ligne par une 
constante : utilisation d'une transformation M. 



Effets de I'arithmetiqu 



Propagation d'errt 



Matrice M 



On divise la premiere ligne de A par 100000 ce qui donne 



0.2000 x10" 4 0.1000 x 10 1 
0.1000 x 10 1 0.1000 x10 1 



0.1000 x 10 1 
0.2000 x 10 1 



en faisant la decomposition LU avec 4 chiffres AVEC permutation on a 

L+U-l 

La solution par LU est "proche" de la solution reelle 



0.1000 x 10 1 0.1000 x 10 1 
0.2000 x10" 4 0.1000 x10 1 



0.1000 x 10 1 
0.1000 x 10 1 



Effets de farithmetique flottante Propagation d'erreurs 



Definition 3.7 



Une mise a I'echelle d'une matrice A consiste a diviser chaque ligne 
de A par le plus grand termes en valeur absolue de la ligne 
correspondante. 



Recette 2 



On peut eviter une partie des erreurs en faisant une mise a I'echelle 
systematique d'une matrice avant la decomposition. 



Attention au determinant 



L'introduction de matrice Maun effet important sur le determinant. Si 
on multiplie une matrice par un scalaire on recupere son determinant 
en divisant le determinant de L par ce scalaire. 



Conditionnementd'une matrice 



On veut etudier le comportement de la solution calculee par rapport a 
la solution exacte pour un systeme lineaire. On va d'abord introduite 
une mesure de la distance dans I'espace R". 



Definition 3.8 * 


Une norme vectorielle est une 


application de R n 


dansR notee 


verifiant 








Q Positivite stricte : 








||x||> 


Vxe 


R n ,x/0 




OVaeK 








\\ax\\ = 


= \a\\\x\\ 


VxeR" 




O Inegalite du triangle 








ll* + y||< 


ll x ll+lly| 


Vx,y g 


R n 



Math-lnfo-ScHS (ENSAM) 



Conditionnementd'une matrice 



Definition 3.9 



La norme euclidienne d'un vecteur x notee ||x|| 2 est definie par 



|X|| 2 = Ixf+xf + .-.+X^ 2 



evidemment la norme euclidienne est une norme vectoriel. 



Definition 3.10 ^ 


II existe d'autre norme sur E". 




» La norme l-\ : 


\\x\\ 






• La norme 1^ 


\\x\ 


oo= max 
/=1 


*/| 



Conditionnementd'une matrice 



Definition 3.11 ^ 


On peut aussi introduire 


une norme 


matricielle. On I' 


a definie comme 


une application A :->• \\A 


verifiant 






Q Positivite stricte : 










\\A\\>0 


VA A + 




O Vael 










\\aA\\= 


a|||/\|| VA 




Q Inegalite du triangle 








ll^ + S||< \\A\\ 


+ \\B\\ VA 


Sde dimensions 


compatibles 


o 








\\AB\\< \\ A \\ 


\B\\ VA Bde dimensions compatibles 


Une norme matricielle DOIT verifier 1-4. 





Conditionnement d'i 



Norme induite 



On peut construire une norme matricielle a partir d'une norme 
vectorielle : 

\\A\\= sup \\Ax\\ 

Une telle norme est dite norme matricielle induite par la norme 
vectorielle. On dit aussi que cette norme est subordonnee a la 
norme vectorielle 



Toute les normes ne sont pas induite 



La norme de Frobenius definie par 



EE1 

est une norme matricielle, s'apparentant a la norme euclidienne mais 
elle n'est induite par aucune norme. 



Math-lnfo-ScHS (ENSAM) 



Conditionnementd'une matrice 



Theoreme 3.5 



Les normes induites par les normes /1 et 1^ sont 

n 

9 \\A\\-\= sup ||/Ax||-i= max Yl IAy| 
||x|h=i 1</<"/=i 



a HA 



sup ||Ax| 

llyll —1 



oo= max £ |Ay| 

1<'<"y=1 



La norme ||>4||i consiste a trouver la colonne dont la somme des 
elements (en valeur absolue) est la plus grande. 

La norme ||>4||oo consiste a trouver la ligne dont la somme des 
elements (en valeur absolue) est la plus grande. 



Conditionnementd'une matrice 



Norme induite par la norme euclidienne 



La norme induite par la norme euclidienne est definie par 



\A\\ 2 = sup ||Axi| 2 = \p{AA') 

|x|| 2 =1 



T\\ 2 



ou p est la plus grande valeur propre (en valeur absolue) de la matrice 
AA T . 

Pourquoi s'interesser au norme matricielle ? On va manipuler des 
vecteurs resultants de produit matrice vecteur. Ce qui nous amene a 
une definition fort importante 



Conditionnementd'une matrice 



Definition 3.12 



On dit que la norme matricielle et la norme vectorielle sont 
compatibles si 

||Ax||< ||>4||||x|| 

pour toutes matrice A et vecteurs x de dimensions coherentes. 

Attention au melange de normes : a gauche norme vectorielle, a droite 
matricielle et vectorielle (respectivement). 



Conditionnementd'une matrice 



Normes compatibles 



On peut montrerque 

• ||A*||i< ||.A||i||x||i la norme 1 des matrices et vecteurs est 
compatible 

• ||^||oo< Halloo IMI00 la norme 00 des matrices et vecteurs est 
compatible 

• ||A>f||2< ||/4||p||x|| 2 la norme Frobenius des matrices est 
compatible avec la norme euclidienne. 

La derniere relation montre qu'il n'y a pas de lien entre norme induite 
et norme compatible. 



Conditionnementd'une matrice 



Soit le systeme Ax = be\ x* une solution obtenue par arithmetique 
flottante. On veut etudier I'erreur par rapport a la solution exacte : 



Si tout va bien ||e||= ||x - x*|| est petit. 

Posons r = b - Ax* r est le residu de notre resolution. 

r = b- Ax* = Ax- Ax* = A(x - x*) = Ae 

On cherche a avoir une borne inferieure et superieure sur ||e||. On va 
se servir de la compatibilite des normes (on choisi une norme 
compatible) 



\A-'r\\< \\A- 



On a aussi 
Cequi nous donne 



l>4e||< \\A\\ 



< II ell < \\A~ 



Conditionnementd'une matrice 



< lle||< 11/4" 



On veut faire apparaitre x dans cette relation (on voudrait Terreiir 
relative") alors on revient au systeme : 



I Axil < IIAIIIIxl 



1 \\A 
— < 



HyII— \\A~' i h\\< II4 -1 llllhll— 
En combinant ces deux inegalites on a 



\h- A \\\\b\ 



< 



1 1 \\A\ 

14-1 llllhll ~~ y ~~ hi 

r\ \\ \\kj\\ -^ *-* 



Conditionnementd'une matrice 



On multiplie la derniere relation par ||e| 



HpII 1411 

<y<uy 



et on a 



Definition 3.13 



4- 1 ||||fo|| " ||x|| " \\b\\ 

A^ \\\\A\\\\r\\ 



IkII ll e ll 

I >A 1 1 1 1 >A — 1 llllhll ~~ llyll — 



Le conditionnement d'une matrice A note condA est defini par 

condA = ||/4" 1 ||||/4|| 



Conditionnementd'une matrice 



Remarque 



» condA depend de la norme matricielle utilisee 
• Puisque 1 < ||/|| parce que 

\\A\\= \\AI\\< \\A 

on a 



et 



1 < ||/||< ||A4" 1 < ||^" 1 ||||^| 
1 < condA < oo 



Theoreme 3.6 



1 \\r\\ llell _ .. Ilrll 
_U — !_!_ < iLJi < CondA— 

condA \\b\\ " \\x\\ ~~ \\b\\ 



Conditionnement d'l 



Concernant le conditionnement 



1 



condA \\b\\ 



< 



\e\ 
\x\ 



< CondA-, 



O le terme du milieu represente I'erreur relative 
Q Si condA est pres de 1 alors 



|e| 
|x| 



\r\\ 
Tb\\ 



est la solution est bonne si ||r|| est petit. 

Q si condA est grand alors I'erreur relative peut etre grande. 

O meme si ||r|| est petit il est possible que I'erreur soit grande 
(manque le conditionnement petit) 

9 Plus condA augmente plus I'algorithme de resolution doit etre 
precis (pour reduire ||r||) 



Math-lnfo-ScHS (ENSAM) 



Conditionnementd'une matrice 



Concernant le conditionnement 



1 



condA \\b\\ 



< 



\e\ 
\x\ 



< CondA] 



Q Meme si condA est petit si la resolution est imprecise on pourra 
avoir des resultats errones (manque ||r|| est petit). 

Q On peut montrer que 



max 



lellllbl 



\r\\\\x\\ 

\e\\\\b\ 



< condA 



O condA depend de la norme matricielle choisie (generalement la 
norme oo) mais ne depend pas de I'algorithme de resolution du 
systeme. C'est une mesure de la sensibilite de la matrice au 
perturbation lors de manipulations matricielles. 



Conditionnement d'i 



Erreurs dans A 



On veut resoudre le systeme Ax = b mais a cause d'erreur de 
troncature ou dans la mesure des elements de A on resout plutot le 
systeme 

(A + E)x = b 

Soit x* solution du systeme perturbe et x la solution du systeme 
original. 

x = /T 1 b = /A" 1 {{A + £>*) = (/ + /A" 1 E)x* = x* + /\" 1 Ex* 

Done si les normes sont compatibles 

||e||= ||x-x*||= ||>!\- 1 £x*||< ||/A- 1 ||||E||||x*|| 

Et on a 

lk-^*ll ,J|£1I 

^ ^ < COndAi--rrr 



Math-lnfo-ScHS (ENSAM) 



Conditionnementd'une matrice 



Erreurs de troncature sur A 



< condA- 



\\x—x*\ 



approxime I'erreur relative 



9 O ressemble a I'erreur relative sur les entrees de la matrice A 

9 Si condA est petit et ||£|| est petit alors I'erreur relative sera 
petite : bonne approximation de x 

9 Si condA est grand on ne peut rien dire 



— ' 



Conditionnement d'l 



Erreurs dans A 



\x — x*\\ ,.\\E\ 

^—r, r— - < COndA-i--rr, 



9 Si E represente I'erreur de troncature sur une machine alors 

\Ejj\ < e m |/\,y| 
par definition de la troncature done 

/ J | t;y| < Cm y j |/\;y| "v4> 1 1 fc 1 1 CXD S Cm||>^||oO 



d'ou 



x - x* 



< tmCondA 



Si le conditionnement est mauvais il faut plus de precision sur les 
nombres 



Math-lnfo-ScHS (ENSAM) 



Conditionnement d'une matrice 



Quoi faire de condA? 



On c'est donne cette mesure qu'en fait-on ? 

9 Exemple 3.16 (3.24) Le conditionnement de A = 270.51 
correspond a une matrice sensible au perturbations. Ne veut pas 
dire que je n'aurai pas la bonne solution mais que j'ai de bonne 
chance d'obtenir une solution erronnees sans un bon algorithme 

» Exemple 3.18 (3.26) condA = 4 ce qui est faible indiquant que si 
tout va bien j'aurai une bonne solution. Rappelons que sans 
permutations la solution est mauvaise ce qui ne contredit pas le 
conditionnement mais indique que mon algorithme est mauvais. 
En incluant les permutations on obtient une bonne solution. 



1. 
1.1 



2 
2 



-10. 
-5.5 



10. 
5. 



cond A = 62 indiquant un mauvais conditionnement. 



Conditionnement d'une matrice 



Quoi faire de condA? 



Est-ce que je peux diminuer le conditionnement en changeant de 

norme ? Non, on a "equivalence des normes" faisant en sorte que 

rien ne changera en prenant d'autre norme. 

Mors quoi? 

» Mauvais conditionnement == appliquer les recettes 1 et 2 
OBLIGATOIREMENT et etre critique des resultats obtenus et 
choisir avec soin I'algorithme de resolution (controler la qualite 
autrement, de maniere ad hoc) 
a Bon conditionnement == les recettes ne sont pas necessaires, on 
peut croire que tout va bien mais on doit encore etre critique des 
resultats obtenus (controler la qualite autrement, de maniere ad 
hoc). 



Conditionnement d'une matrice Retour 



9 Decomposition LU 

9 Structures particulieres menant a des algo simplifies et plus 
efficaces : 

o A symetrique, definie poisitive : Choleski 
• A tridiagonale 

9 Effets de la representation virgule flottante 

« Recette 1 : permutation des lignes apres chaque construction de 

colonne de L 
» Recette 2 : mise a I'echelle des lignes de la matrice 
» Conditionnement : un indicateur de la sensibilite du systeme lors de 

la resolution 
<» CondA + r = (Ax - b) donne un verdict positif si les deux sont petits 
<» CondA seul = rien 
a r — Ax - b seul = rien 



Methodes iteratives pour la solution d'un systeme lineaire Principe des methodes iteratives 



Exemple introductif 

Supposons que nous cherchions a resoudre le systeme suivant : 



2*1 


+4x 2 


-2x 3 


-6x4 


= -6 


1*1 


+3x 2 




-6X4 


= 


3x, 


-x z 


+*3 


-6X4 


= 8 




-x 2 


+2x 3 


-6x4 


= 6 



et que nous connaissions un vecteur x*i , x* 2 , x* 3 , x* 4 proche de la 

solution. 

-6 - 4x* 2 + 2x* 3 + 6x* 4 
x t 

•*3 

X 4 



0- 


-1X*i 


2 
-0x* 3 + 6x% 


8- 


-Ax*, 


3 « « 

+ X* 2 +6X* 4 


6- 


-Ax*, 


1 
+ X* 2 -2X*3 



La politique des « petits » pas. 



on cherche a creer une suite de vecteurs {x^ k \k e N} qui converge 
vers x solution de Ax = b 



Math-lnfo-ScHS (ENSAM) 



(*) 



* Jacobix) ',k e IN 



(*) 



9 Gauss Seidel x££, k e N 
9 Relaxation x^\ k e N 



a chaque iteration 



le vecteur d'erreurs : e^ = xW - x doit diminuer 



J 



comment choisir? 

• il en existe plein d'autres... 

9 la question centrale : quand convergent'elles ? 

• seconde question : a quelle vitesse 



parmi celles qui convergent on prend celle qui converge le plus vite ! 



Methodes iteratives pour la solution d'un systeme lineaire La methode de Jacobi 



/-1 



Ax = b 



Y^ AijXj + AjiXj + Y A ij x i = b i / = 1 , n 

y=1 y=/+l 



Fonction x «- Jacobi(/\,b,x) 



Tant que (on a pas converge) faire 
Pour /'= 1, a? faire 

"i ~ zJy=i Ay*/ ~~ zJy=/+i Ay*y 



Fait 



Fin Pour 

Pour /= 1, a? faire 

| Xi = V\ 
Fin Pour 



An 



wm 



Exemple : methode de Jacobi 



3 

1 




avec b 



etx 



X1 



*y 2 = ^ 



• ys 



i 

3 

1+1/3-1/2 



_3_ 
12 



_5_ 
24 



*2 




norm(/4*x -b) = 1.7321 
norm(/A * X1 - b) = 0.4488 
norm(/A * x 2 - b) = 0.1718 



[x 2 A\b] 



'0.2500 0.294 f 

0.3333 0.3529 

,0.2083 0.2353; 



Les iterations de Gauss Seidel 



Fonction x <- Gauss_Seidel(/4,b,x) 

Tant que (on a pas converge) faire 
Pour /= 1, a? faire 

"i ~ Zj/=1 Ay*/' ~ 2^/=/+1 Ay*/' 



Fait 



X, - - 

Fin Pour 



plus « simple », mais non parallelisable 



Les iterations de relaxation 



Fonction x «- Relaxation(/4,b,x,w) 

Tant que (on a pas converge) faire 
Pour /'= 1, a? faire 

x, = ,, i b '-^=!^-^/+i^ x A + (1 _ u)Xl 



Fait 



A„ 



Fin Pour 



pour uj = 1 c'est Gauss Seidel 



Quand s'arreter ? 



^new v old II ^ -. 



» ||x n 

» un nombre maximal d'iterations est atteint 



Fonction y «- Gauss_Seidel(/\, b,x,£, n_ite_max) 
Pour / = 1, n faire 

| y, = *i xi = y, + 2 * £ 

Fin Pour 

n_ite = 

Tant que (||x - y|| > e ET n_ite < n_ite_max) faire 
Pour /' = 1 , n faire 

| Xi = 7/ 
Fin Pour 
Pour / = 1 , n faire 

M-1 



y,= 



"i z2j=i Af/y zJy=;+i Ay*/ 



Fin Pour 

njte = njte + 1 



Fait 



• x 

• y 



Formulation matricielle des iterations 



Considerons maintenant la suite de vecteurs xM construit a I'aide des 
matrices carrees M et N, ainsi que d'un vecteurb : 

MyS k+ ^ = AfeW +b 

On peut reecrire ces iterations a I'aide de la matrice carree C = /W 1 A/ 
et du vecteur d = M~ 1 b : 

X (*+1) = M -i N X W + w -i b 
C x( fc ) + d 

pour un vecteur x^ ) donne. 



□ rS 



Decomposition d'une matrice 



Revenons au systeme Ax = b. Nous pouvons decomposer la matrice A 





A 


= 


D 


+ 


u + 


L 


fa u 


a 12 


a is\ 


fa u \ 




(0 3i2 3i 3 \ 


/ 0> 


321 


a 22 


a 23 J = 


3 22 


+ 


3 23 + 


3 2 1 


\3 3 i 


a 32 


a 33/ 


V 3 33 / 
diagonale 




\0 / 

triangulaire sup 


\ a 31 332 0y 

triangulaire inf 




Ax 




Dx 




+ 


(L+U)x 


a 11 x 1 + 312X2 + 3i3X 3 \ 


/3iiXi +0 




+0 \ / 


+812X2 +813X3 


321 X1 + 822X2 + 823X3 


= +3 22 X 2 




+0 + 3 2 iXi 


+0 +3 23 X 3 


a 31 x 1 + a 32 x 2 + a 33 x 3/ 


V +0 


+ a 33 x 3/ \ a 31 x 1 


+ 3 32 X 2 +0 



Ax = Dx + Lx + Ux 



Decomposition d'une matrice 



Revenons au systeme Ax = b. Nous pouvons decomposer la matrice A 



De maniere generale : A = D + L + U avec : 



D = diag(a u ,a22,... 


Snn) 








/ ... 
a 21 




\ 




L = 


^31 a 32 
\3n1 a„2 


&n,n— 1 







U = 


/0 312 313 




a n-2,n- 


am \ 

a 2n 

1 a n-2,n 




VO 









/ 



□ rS 



Les iterations de Jacobi s'ecrivent alors 



Le point fixe 



Ax =b 

(D + L+U)x =b 

Dx + (L+U)x =b 

Dx 



(L+U)x + b 



Les iterations de Jacobi 

Dx (/c+1) 

x (*+ 1 ) 



-(U+L)xW+b 
-D-\U + L)xW + D~\ 



□ rS 



Les iterations de Jacobi s'ecrivent alors 



Fonction x 4- Jaoobi(A b,x) 

Tant que (on a pas converge) faire 
Pour / = 1 , n faire 



bi-J2 A a x i- } 
y'=i y'=/+i 



Fin Pour 

Pour /' = 1 , n faire 

I ><i = y, 

Fin Pour 



Fait 



Fonction x •<— Jacobi(/A, b.x) 


Tant 


que (on 


a pas converge) 


faire 






I 


c=D\(fc- 


-(L+U)x) 


Fait 







,(*+1) 



■'y 



xW + 



-v 



-D-1(U+/.) 



D- 1 b=D\fc 



Gauss Seidel et la relaxation 



les iterations de Gauss Seidel s'ecrivent 

(D + Z.)x(* +1 ) =-L/xW+b 

x (/c+i) =(D + /.)- 1 (-L/)xW + (D + /.)- 1 b 

et les iterations de relaxation s'ecrivent 

(D + w^x^ 1 ) = ((1 -^D-wl^xM+cjb 

x (*+1) = (( D + w/ _))-1 ((-i _ W ) D _ wL/ ) X W + 



c r 



{(D + uL))-\ojb) 



Pour w = 1 , on retrouve bien I'equivalence des formulations de Gauss 
Seidel et de la relaxation 



□ [5 



Condition suffisante de convergence 



X (*+1) = CfcW + d 
Si x est solution du probleme, on a x = Cx + d et done 

x< k + 1 )-3? =C(xW-x) 

= C 2 (x(' ( - 1 )-x) 

= C* +1 (x(°'-x) 

Si maintenant on raisonne en norme : 
erreur a I'etape /<+i 



II* 



(*+i) 



= ||C /( + 1 (x( )-x)|| 

< IIC^IIIIxW-xll < ||C|| /<+1 llx^-xl 



-►o ? erreur initiale 
ce faisant, nous avons demontre le resultat suivant : 



Theorem (convergence d'une suite vectorielle) 



S'il existe une norme matricielle telle que \\ C\\ < 1 alors la suite xM 
converge vers x = (I - C)~ 1 d 



oaj 



Normes matricielles de type vectorielle 



Norme vectorielle (norme de Frobenius) 



|C||| = ££c 



Ce type de norme n'est pas toujours pertinent pour ce que nous 
cherchons a faire. En effet, certaines normes vectorielles ne sont pas 
sous multiplicatives (la propriete \\AB\\ < ||/\||||6|| n'est pas toujours 
verifie) ce qui est facheux pour notre analyse. 
Par exemple pour la norme vectorielle \\A\\ A = max/j \Ay\ : 



A = B 



1 1 
1 1 



\AB\\ 



2i ||^|| A ||fi||A = 1 



Normes matricielles d'operateur 



Pour conserver cette propriete, il faut utiliser une Norme d'operateur 
definie par : 

||C|| = sup 

x^O ll x ll 



Par construction on a bien : 

\\ABx\\ < 11/411116x11 < ||/A||||ei|||x| 



=► \\AB\\ < \\A\\\\B\\ 



Exemples (theoremes) : 



|C||i 

\r\\ 



max 

J 



£ic, 



max 

/ 



£i<* 



|C|| 2 = 
» pour C symetrique : max |A, 



ou fn est valeur propre de C T C. 
A; = JJij est val. propre de C. 



demonstrations : de maniere generale nous allons montrer les deux 
inegalites suivantes 

<P(C) < \\C\\ < <j>{C) 

On commence par trouver une borne sur cette norme : 

|| C | h =supi3 

x^O ll x ll 
= sup '-— '- 

Xy^O |W| 

. E/Ey|Cy||x y | 
< sup -j— 

x^O ll x ll 

<supE(EI^I)S1 

x^O y ,■ ||X|| 

<sup£(max£|q|)^ 
X^O y y ,■ 11*11 

<sup(max^|C^|)^M 

X^O J ; y l|X|| 

<max^|Cj| car m <1 

en choisissant le vecteur e y = (0, . . . , 0, 1 , 0, . . . , 0) T , on atteint la borne 



Le cas de la norme 2 



on passe par le systeme des vecteurs propres : C T C\j = ^,v, 

x = X>' «- x = P£ ; || X || 2 =^ T P T PC = ^^ 

/ / 

IPyII 2 

l|C||l = sup^— "]- 

x T C T Cx 
= sup 

x^o IWr 

= sup — ^p - 

x^o INI 

= sup xT( p,g WY,) 

xyo IWr 

= sup ' ' < max/i, 

x^o INI ' 

La encore, en choisissant le vecteur propre associe, on atteint la 

borne : 

HCv,|| 2 _ T _ 

n no = V/C CV/ = /i/ 

M 



Illustration 2d. 



\A\\ = sup 



1/4x1 



sup \\Ax\ 



x/0 ll x ll ||x||=1 



Definition 



rayon spectral 



p(A) = max I A,- 1 

1</<n 




Proprietes des normes (Op) ^ Proprietes des normes 



• ||i4|| > Osi A / 

9 ||>A|| =0o A = 

9 \\kA\\ = \k\\\A\\ 

9 \\A + B\\ < \\A\\ + ||B|| 

o \\AB\\ < \\A\\\\B\\ 



9 \\A\U =maxyE/IAyl 
» \\A\loo = max,-X;ylAyl 
9 Asym. : \\A\\ 2 = p{A) 
9 A sym. : p(A) < \\A\\ 

• \\A\\ 2 < VIAhl^lloo 



Resumons nous (normes matricielles?) 



Q on cherche a resoudre Ax = b 

Q on utilise une methode iterative : x^ +1 ) = Cx^ +d 

O I'erreur... e W = C k e^ 



O ... peut etre controlee : 



;W|| < II CI 



,(0)| 



|Cx|| 



O si I'on trouve une norme telle que ||C|| < 1 la suite converge 
Q definir une norme telle que 
Q definition : 

O theoremes (bien pratiques...) : 

» llClh = max J2\Cij\ 



\C k \\ < \\C\ 
= sup x ^o ■ 



|C||i 

IdloC 



i 

max^|C^ 
i 



* \\C\\z = max- 
Q on teste les differentes normes pour avoir ||C|| < 1 

A quelles conditions sur A, la methode iterative converge ? 



J 



Le cas des matrices a diagonale dominante 



Definition 




Theorem 



Si A est a diagonale dominante alors les methodes de Jacobi et de 
Gauss Seidel convergent 



Demonstration : 



9 Jacobi 



J J \\oo 



\D~\M + N)l 



max . , 

/ IA 



oElA, 



< 1 



i+i 



9 Gauss Seidel : on montre que si Gauss Seidel diverge, alors A 
n'est pas a diagonale dominante : 

1 < max,- 4j Y^'jZ] lAyl + S"=/+i IAy|- Les details de la preuve ne 
sont pas evidents... 



Exemple : la matrice est-elle a diagonale dominante ? 



'3 -1 r 

1 -2 o 

.1 1 4, 



( 2 


-1 





°\ 




/-9 3 1 


-1 




2 

-1 


-1 

2 




-1 


;^ 3 = 


1 2 
-1 1 8 


Vo 





-1 


2; 




\-2 1 4 



Remarques : 

9 toute matrice a diagonale dominante est reguliere 

9 on retrouve ce type de matrice dans les methodes a elements 
finis (c/meca/EP) 



□ rS 



Le cas des matrices positives 



Definition (rappel) 



Une matrice carree A est definie positive si elle est symetrique et si 

Vx / G R" x T Ax > 



Theorem 



Si A est strictement symetrique definie positive alors la methode de 
Gauss Seidel converge et la methode de la relaxation pour < w < 2 



□ [5 



Etude des perturbations (ou des erreurs) 



/10 7 8 7\ 




/32\ 


7 5 6 5 


■ b — 


23 


8 6 10 9 


i ** — 


33 


\7 5 9 10/ 




V3iy 



admet comme solution x 



/ 1 \ 
1 

1 

vv 



b + S b 



/32,1\ 
22,9 
33,1 

\30,9/ 






admet comme solution x + 5 X 






comment mesurer cet effet ? 



Etude des perturbations (ou des erreurs) 



/10 7 8 7\ /32^ 

7 5 6 5 23 

8 6 10 9 ; 33 
\7 5 9 10/ \31/ 



admet comme solution x 



/ 1 \ 
1 

1 

W 



b + S b 



/32,1\ 
22,9 
33,1 

\30,9/ 



admet comme solution x + 5 X 



( 9 ' 2 \ 
-12,6 

4,5 

V 1,1 I 



Pour la norme infinie : 

ll^llc 



0,003 



ll&fl 



13,6 



C'est la nature de la matrice A qui est en cause : 
comment mesurer cet effet ? 



Conditionnement d'un systeme lineaire 



Ax = b 

9 ||b|| = \\Ax\\ < \\A\\ ||x|| 
9 \\5 X \\ = \\A~'5 b \\ < \\A-'\ 

IIM 



A(x + 6 X ) = b + 5 b 



1 \\A\ 
x - b 



IN 



IM <II*~ 1 II \\6 b \ 



< n/i-1 1| WAW^l 
xii \\b\\ 



cond(A) : conditionnement de A 




9 cond(/4) proche de 1 : stable 
9 cond(A) grand devant 1 : instable 



□ rS 



Illustration en 2 d 



Ax 



a u x^ +a 12 x 2 = bi 
321*1 +a 22 *2 = b 2 




cond(/\) proche de 1 : stable cond(/4) grand devant 1 : instable 



Methodes iteratives pour la solution d'un systeme lineaire 



Introduction 

9 Introduction 

Elimination de Gauss 

» Generates 

9 Methode des substitutions successives 

» Matrices diagonales 

9 Matrices triangulaires 

9 Algorithmes de resolution pour les matrices triangulaires 

9 Unicite dans le cas des matrices triangulaires 

9 Exemple de resolution 

9 Comment exploiter les algorithmes des matrices triangulaires? 

9 Transformation des matrices 

9 Operations elementaires sur les lignes 

9 Operations sur les lignes : matrice M 

9 Operations sur les lignes : matrice P 

9 Operations sur les lignes : matrice T 



Ho O,: 



Math-lnfo-ScHS (ENSAM) 



Conclusion 



9 3 methodes iteratives (Jacobi, Gauss Seidel, relaxation) 

9 chaque iteration < 0(n 2 ) : produit matrice vecteur 

• encore mieux si A est sparse 

» peut s'appliquer a des matrices par blocs 

• convergent a certaines conditions sur A 

9 pourcomprendre : il faut connaitre la notion de norme matricielle 

• I'importance des valeurs propres 

9 il existe des methodes plus sophistiques 

» a chaque iteration on se rapproche de la cible... 



□ si 



Differents types de matrices 



» quelconque : n lignes p colonnes 
• petite, moyenne, grande 
a pleine, creuse, par blocks 
a complexe, reelle, binaire 



QR 



» carree : n lignes n colonnes 

» quelconque (carree) 

a symetrique 

a symetriques positive 

a diagonale dominante 

a triangulaire 

a diagonale 



Differents types de matrices 



» quelconque : n lignes p colonnes 
• petite, moyenne, grande 
» pleine, creuse, par blocks 
a complexe, reelle, binaire 



QR 



» carree : n lignes n colonnes 

9 quelconque (carree) LU 

a symetrique LDL 

a symetriques positive Cholesky : LL', GS 

» diagonale dominante Jacobi 

9 triangulaire trisup et trijnf 

9 diagonale